/*This file will calculate the test score gap under a fixed rank
and distribution as chosenin lines 18 and 20*/



clear all

/*Choose data file here*/
cd "C:\Users\bond10\Desktop\RA\Test Scores\Paper\replication files"

use ECLSbondlang.dta

set more off

mata

/*what rank ordering*/
st_view(k1=.,.,("readfallk weights black white"))
/*what distribution*/
st_view(t1=.,.,("readfallk weights black white"))

k1 = sort(k1,(1))
t1 = sort(t1,(1))
rows(t1)

currenttarget = 1
sumblack = 0
sumwhite = 0
currenttargetweight = t1[currenttarget,2]
sumweight=0

for (i=1; i<=rows(k1); i++)   {
      currenthostweight = k1[i,2]
	  while (currenthostweight>0) {
	  if (k1[i,3]==1) {
	    if (currenttargetweight<=currenthostweight) sumblack = sumblack + t1[currenttarget,1]*currenttargetweight
		else sumblack = sumblack + t1[currenttarget,1]*currenthostweight
		if (currenttargetweight<=currenthostweight) sumweight = sumweight + currenttargetweight
		else sumweight = sumweight + currenthostweight
		              }
	  if (k1[i,4]==1) {
	  if (currenttargetweight<=currenthostweight) sumwhite = sumwhite + t1[currenttarget,1]*currenttargetweight
	  else sumwhite = sumwhite + t1[currenttarget,1]*currenthostweight
	  		if (currenttargetweight<=currenthostweight) sumweight = sumweight + currenttargetweight
		else sumweight = sumweight + currenthostweight
                     }
	  usedhostweight = currenthostweight
	  usedtargetweight = currenttargetweight
	  currenttargetweight = currenttargetweight - usedhostweight
	  currenthostweight = currenthostweight - usedtargetweight
	  if (currenttargetweight<=0) {
	  currenttarget = currenttarget + 1
	  currenttargetweight = t1[currenttarget, 2]
	  }
	  }
        }
blackweight = k1[.,2]'*k1[.,3]
whiteweight = k1[.,2]'*(k1[.,4])

sumweight

mb = sumblack/blackweight
mb
mw = sumwhite/whiteweight
mw

weights = t1[.,2]
sum(weights)
mu = (weights't1[.,1])/sum(weights)
mu
demean = t1[.,1]:-mu
top = weights'(demean:^2)
sdev = sqrt(sum(weights)*top/(sum(weights)^2-sum(weights:^2)))
sdev

diff = (mw-mb)/sdev
diff

end
